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^ ' Abstract 
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k> , We present a multi-grid algorithm in order to solve numerically the thermody- 

namic Bethe ansatz equations. We specifically adapt the program 
data of the conformal field theory reached in the ultraviolet limit. 
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PROGRAM SUMMARY 

Title of program: TBA 

Catalogue number: 

Program available from: CPC Program Library, Queen's University of Belfast, N. Ireland 

Computer: IBM RS6000/560 

Operating system: AIX 3.2 

Programming language used: FORTRAN 77 

Number of lines in program: 1287 

Key words: two-dimensional systems away from criticality, conformal field theory, ther- 
modynamic Bethe ansatz, multi-grid methods. 

Nature of the physical problem: resolve the thermodynamic Bethe ansatz equations cor- 
responding to a factorized scattering theory and extract the data of the underlying 
conformal field theory reached in the ultraviolet limit. 

Method of solution: non linear multi-grid method with full adaptive scheme. 

Typical running time: strongly dependent on the dimension of the system of coupled 
integral equations, see Table 1. 
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LONG WRITE-UP 

1 Introduction 

Massive relativistic field theories can be described on-shell by their scattering matrix. 
This approach is specially fruitful in two dimensions, where there exists a large class of 
models which are integrable, and their S-matrix can be computed exactly, being fac- 
torizable [lj. Unfortunately there is no general direct method in order to compute the 
S'-matrix of a theory, but usually it is conjectured from general axioms and the underly- 
ing symmetries of the corresponding Hamiltonian. 

The thermodynamic Bethe ansatz (TBA) was developed in order to provide a means to 
link a conjectured scattering theory with the underlying field theory B. It describes the 
finite temperature effects of the factorized relativistic field theory, using the S matrix 
as an input. If one studies the high temperature limit of the TBA equations, one can 
identify the conformal field theory (CFT) which governs the ultraviolet behaviour of the 
underlying field theory. One should though note, that it is not guaranteed that every 
consistent S'-matrix describes the scattering in some field theoretical model! Therefore 
the axiomatic bootstrap approach is only of limited value if not linked to field theory by 
some means, wherefrom the TBA is one of the most powerful ones. 
Given the scattering data one can in most cases extract analytically the central charge of 
the CFT reached in the conformal limit, and in some cases the dimension of the perturb- 
ing operator, if the symmetry of the problem is known. Numerical calculations on the 
other hand can solve the TBA equations and therefore extract any measurable quantity. 
In 0, U the TBA equations were resolved by an iterative method. We propose here a 
multi-grid algorithm, which is considerable faster, an important fact if many particles are 
involved. The heart of the program is the resolution of the coupled integral equations. 
Around this core we have designed some utility-programs, in order to make the tool easier 
to use. We specialize our application to the case of diagonal S-matrices, see e.g. |2|, |3|. H]. 
As physical quantities we extract the central charge, the dimension of the perturbing 
field and the perturbation expansion. Note though, that one can easily add subroutines 



calculating other quantities. 

2 The TBA Equations 

Consider an integrable massive scattering theory on a cylinder. This implies factorized 
scattering, and so one can assume that the wave function of the particles is well described 
by a free wave function in the intermediate region of two scattering. Take the ansatz 

V(xi...z n )=e i I>i*i£A(P)e(zp) ; 

p 

A(P) are coefficients of the momenta whose ordering is specified by 

e(i f 1 if x pi <...<x p „ 

I otherwise 

Let the permutation P differ from P' by the exchange of the indices k and j. Then 

A(P') = S kj ((3 k -(3,)A(P) ■ (1) 

We impose antiperiodic boundary conditions for our wave functions, which provides that 
two particles cannot have equal momenta, leading to the condition 

A(k,p 2 ,..., Pn ) = -e i ^ L A(p 2 ,..., Pn ,k) , (2) 

L being the length of the strip on which we consider the theory, comparing (0) and (0), 
one realizes 



iLmkSinhPk U S ki(Pk-P J ) = -l for fc = l,2,...n. (3) 



e 



We introduce the phase 5kj(Pk ~ Pj) — ~ i hi Skj(Pk ~ Pj)- in terms of these the equation 
become 

Lm k sinh (3 k + ^ hj{Pk ~ Pj) = 2vrn fe for k = 1, 2, . . . , n , (4) 

rik being some integers. These coupled transcendental equations for the rapidities are 
called the Bethe ansatz equations. One tries to solve these equations in the thermody- 
namic limit introducing densities of rapidities for each particle species and transferring 



the equations into integral equations. That is, let pf {(3) = xg, where we assume that 
there are n particles in the small interval A/3, be the particle density and p^ a \(3) = ^% 
be the level density corresponding to the particle a, then (ffl) becomes 



m a Lcosh(3 + J2 f°° Vab{f3-P')Pi\f3'W = 2tt P ^ . (5) 

6=1 J -°° 

In order to compute the ground state energy one needs to minimize the free energy 

RLf(p,p 1 ) = RH B (pi) + S(p,pi) , (6) 



where Hb = So rn a I cosh f3p\ d(3 and 5* denotes the entropy. The extremum condition 
for a fermionic system^ takes the form 

dp 

J fab{P ~ P ) iogv f e -*■') 

6=1' 



n r°° , , dff 

~M a coshf3 + e a (f3) = Y, <P<*(P ~ ^)log(l + e~^)^- , (7) 



where we introduced the so-called pseudo-density e 6a = , * (a) , the scaling length 

py a >— p\ 

r = Rmi and the rescaled masses M a = — ; mi is the lightest particle mass. These 
coupled integral equations are called the TBA equations. The extremal free energy 
depends only on the ratios ^j and is given by 

f( r ) = -7T E M W cosh/?log(l + e- £ ^))d/3 . (8) 

2vr a=1 J-oo 

One can extract several physical quantities from the solution of the TBA-equations 

([fj. |], |3|]). Since very little is known about non-critical systems, one tries to examine the 

equations in the ultraviolet limit, which corresponds to r — ► 0, where the underlying 

field theory should become a CFT. The central charge is related to the vacuum bulk 

energy, and is given by 

Qy n />00 

c(r) = — VM a / cosh/?log(l + e~ ea(/3) )c^ . (9) 

71 a=i J -°° 

Having calculated the central charge one would like to extract the conformal dimension 
of the perturbing operator. For small r, one expects that /(r) reproduces the behaviour 
predicted by conformal perturbation theory, which in terms of c(r) reads as 

Of OO 

c ( r ) =c _^V + £ /fcr y* 9 (10) 



1 We use the fermionic TBA equations since in diagonal scattering up to now they turned out to be 
the relevant ones, see e.g.H for the general theory 



with possible logarithmic corrections. 

The exponent y is related to the perturbing field by y = 2(1 — A) if the theory is 
unitary and by y = 4(1 — A) if it is non-unitary The coefficients are related to correla- 
tion functions of the CFT |2|, [|, and even if one cannot read them off directly, this is an 
ultimate important check of the theory. 

Note that the application chosen is not a limitation of the use of the program. Also 
non-diagonal S-matrices (see ||) can be treated, since once one has diagonalized the 
transfer-matrix also in that case the numerical problem reduces to solving fl7|). Further 
quantities to measure can simply be added, and also one can study any range of r, being 
a parameter in the input-data. 



3 Description of the Solution Method 

Multi-grid (MG) schemes are known to be the most efficient methods for solving elliptic 
boundary value problems. Actually, the underlying idea of treating the different char- 
acteristic length scales of the problem on different grids, applies successfully also to the 
numerical resolution of various other problems, as the resolution of integral equations 
H, |n[ . The system of non linear Fredholm integral equations (0) has been solved using 
iterative methods [§, [3|. Even if these methods provide a satisfactory solution in terms 
of accuracy, the number of iterations and corresponding computer process (CPU) time 
required to reach a specified precision can become excessively large as the number of grid 
points N increases. Typically a simple one level relaxation would require 0(N 2 log N) 
operations. With a multi-grid solution technique the computing time for integral equa- 
tions is reduced to 0(N 2 ) |TT[, and in particular cases to O(NlogN) ||, thus justifying 
the extra effort in programming. 

Now we define our numerical problem and we explain how the multi-grid scheme 
works for solving it. In discretising the TBA equations (0), we use the trapezoidal rule 



flQfl on a grid with mesh size h so that our system yields 

h n 

e a (/3) =rM a cosh /? + —£ £ w (f3')pabW ~ P')log(l + e^)) , (n) 

llX b=l/3'eCl h 

a = 1,2, ... ,n, (3 G Qh, where Qh is the set of grid points with grid spacing h. The 
weights are w(f3) = 1 unless on the boundary where w(j3) = 1/2. Now let us introduce a 
sequence of grids with mesh sizes hi > h% > ... > Km-, so that h^\ = 2h^. The system 
fllTl) with discretisation parameter hi will be denoted as 

< = Kf*(4) + fi » o = l,2,. ..,n , (12) 

where a summation over b is intended and where 

Ki b (e)(/3) = ^- £ «;(^)^(/3 - /3') log(l + e-^')) . (13) 

/7r /3'en fe , 



Following |TT[ we have applied one (Gauss-Seidel) iteration to ([12]), and obtained the 
approximated solutions e l a , a — 1, 2, . . . , n. We then transfer them onto the next coarser 
grid, e 1 ^ 1 = If 1 ^ where l[~ x is a restriction operator. The coarse grid equations become 

& 1 = K^fr 1 ) + fi~ X > a=l,2,...,n , (14) 

where 

Ja =H Ja + e a ~ K ab i e b ) ~ H ( 6 a ~ K ab{ e b)) > ( 15 ) 

and with J^ _1 another fine-to-coarse grid transfer operator not necessarily equal to I^- 
Having obtained the solution of the coarse grid equation e l ~ x the difference e l ~ x — e^ 1 is 
the coarse-grid (CG) correction to the fine-grid solution 

e a <- e a + ^-l( e a ~ e a ) , (16) 

a = 1, 2, . . . , n, and j|_ 1 is a coarse-to-fine grid interpolation operator. Finally we perform 
one relaxation at level £, in order to smoothen errors coming from the interpolation 
procedure. To solve the system of equations (0) we employ a coarse-grid correction 
recursively, i.e. equation ([14]) is itself solved by iteration sweeps combined with a further 
CG correction. 



4 Numerical Performance 

The algorithm used to solve equations ([12]) is a non-linear multi-grid (NMGM) method 
( ||1 1|| ) with full adaptive scheme (FAS) (||). The program can perform different MG 
variants: V or W cycle and (non linear) nested iteration, depending on the value of the 
parameter chosen. In order to reduce the number of parameters as input we have fixed 
most of them, that is the parameter regarding the multi-grid cycle (in order to optimize 
it), leaving as an input only the physical ones. 



Figure 1: Evolution of residual error norm with CPU time for a 1-particle system at 
r — 0.1 for different HX: solid line for MG, dashed line for iteration only. 

For any scaling length r we use an initial approximation which behaves like rM a cosh/9, 
wherefrom the program determines the numerical boundary at which the kernels vanish 
and verifies that the conditions for the existence of (at least) one solution given by the 



Schauder's fixed point theorem are satisfied [12). Having determined the size of the 
numerical domain for a given r the number of levels M is set such that the finest level 
has a mesh-size hu of order of HX, which is one of the input parameters. 



no. equations 


CPU time (sees) 
Relax Multi-Grid 


1 


4 


3 


2 


34 


22 


3 


508 


331 


4 


1230 


712 


5 


2530 


1320 



Table 1: A comparison of CPU time required to reach a particular value of the norm 
(19), for r = 0.1, HX = 0.1. 

We compare the performance of the MG and of a (Gauss-Seidel) iterative scheme 
in terms of CPU time in Figure 1, there the different initial residual error for MG and 
iterative scheme is due to the set up of the initial approximated solution in the MG cycle, 
that is a non-linear nested iteration which uses a MG cycle itself (see ]11[). As a norm 
for the residuals 



r (/3) = (e a - K ab {e b ) - f a )((3), (3 e Q hi 



we define the norm 



t || M = max / J2 T M' d 

Ka<n \ „ '—' 



(17) 



;is) 



In order to outline how the multi-grid algorithm becomes important as the number of 
particles increases we give in Table 1 the CPU time required by the two methods to solve 
the discretized problem to a value of the residual norm 



r||M<l-10 



-ii 



(19) 



5 Structure and Use of the Program 

As already mentioned in the introduction the program consists of two parts: the core, 
which resolves the TBA-equations (0) and the periphery, which on the one hand con- 



structs the kernel, and the initial solution , and on the other hand extracts from the 
solution the central charge, the dimension y and the coefficients f\ of the perturbation 
expansion fllCf). 

We specifically designed the program for diagonal scattering theories, that is we are 
concerned with scalar S matrices of the form 

Sab = Y[f(aab) , i = l,-..,n ab , a,b=l,...,n (20) 

i 

with 

sinh |(/3 — ma) 
n being the number of particles in the theory and n a b is the number of factors f x ap- 
pearing in the S'-matrix S a b (for a recent review on this subject and many examples, see 
0). The set of the numbers a, and the masses of the theory are sufficient to resolve the 
TBA-equations. 

As input-data we have therefore three data files: TBA.DAT containing general infor- 
mation about the range of r, the mesh-size and the number of particles, and the file 
ALPHA.DAT containing the values of a and MASS. DAT containing the values of the 
mass of the particles. From the input the program then constructs the kernel, and the ini- 
tial approximation, which are then used by the multi-grid algorithm. The routine returns 
the solutions encoded in a matrix array Q which is stored in the output file SOL. DAT. 
The physical parameters are then calculated in the appropriate following subroutines, see 
the flow-diagram. 

We discuss here specific structure of the algorithm in terms of a simple example. 
Consider the S'-matrix 



Sii = /a/a, S12 = £21 = /1/2/3/4, S22 = fif±(fifi) , (21) 



\2 

55 5555 5555 



being known to describe the minimal model M2J, (i-e. c e // = |) perturbed by the field 



7 ' 



with dimension A = — I. The masses of the two particles are 



Mi = 1 , M 2 = 2cos(-). (22) 

5 



Then the input-file ALPHA.DAT should read as 



2.0 5.0 

3.0 5.0 

-2.0 1.0 

1.0 5.0 

4.0 5.0 

2.0 5.0 

3.0 5.0 

-2.0 1.0 

2.0 5.0 

3.0 5.0 

2.0 5.0 

3.0 5.0 

1.0 5.0 

4.0 5.0 

-2.0 1.0 
Every factor f a is characterized by two numbers which are the nominator and denomina- 
tor of a. As a field separator between values a belonging to different ^-matrix elements 
one gives a number smaller than —1. Since the S- matrix is symmetric the program reads 
only the elements S\^, Sip, ■ ■ ■ , Si :Jl ; 5*2,2; • • • , S2, n ', ■ ■ ■', Sn-i,n-i, Sn-i,n] S rhn . 

The file MASS. DAT contains two values: 
1.0 

1.618033988749895 
Finally the file TBA.DAT contains in our example the following data 

15,2 11,12 

1.0d-14,1.0d-2 ZER0,HX 

0,1 NREL,IWRITE 

corresponding to 

30,0.01,0.01 MAX, STEP, R0 

20.0,7.0,0,0,4 YN,YD,NY,NCEX,MFIT 

4., 7. CEXN,CEXD 

These parameters are: 

11,12: corresponding lengths of the files ALPHA.DAT and MASS. DAT (12 coincides with 



the number of equations of our system, NSYS in the program); 

ZERO: the order of the value of the residual norm to be reached; HX the finest grid size; 

NREL: if NREL= multi-grid algorithm, if NREL= 1 only relaxation is done 

MAX: the number of different radii r to be used; RO is the smallest radius, and STEP give 
the step-size of the radius in the cycle; 

YN, YD: are the numerator and the denominator of the exact exponent if available, otherwise 
dummy numbers; if NY= the exact exponent is used in the fitting procedure, if 
NY= 1 the estimated y is used; if MFIT> the program calculates the first MFIT 
coefficients fc (in this case must be MAX > MFIT + 5). Finally if NCEX= 0, the 
exact charge (cexact=CEXN/CEXD) is used for the fitting procedure. 

Some comments: The CPU time is mainly determined from the parameters NSYS, 
being the dimension of the system, MAX and HX. MAX is chosen corresponding to 
the physical parameters one wants to compute: in order to get a sensible result for the 
dimension y it is enough to use MAX ~ 5, whereas if one wants to calculate the coefficients 
fi a larger number is required; for example with MAX> 10 one can get f\ up to 0(1CT 5 ); 



the more data is used the better the fit-procedure |L3[ works and more coefficients can 
be obtained. HX determines the error of the integration routines calculating c(r) . In any 
case HX=0.01 should be sufficient. 

The program produces four files of output. If IWRITE=1 only the physical informa- 
tions, that is c, y and fi are stored in the file OUTPUT.DAT (see test-run output). When 
IWRITE=0 the program generates, in addition, the file RES. DAT which containes tech- 
nical data and the evolution of the value of the residuals during the iterations. Together 
with RES. DAT other two files are written: TIME. DAT which containes the logarithm of 
the residual norm and the actual CPU time used, and SOL. DAT (where the solutions for 
any r are stored). 

6 Test-Run Output 

The following file was produced using the input files described in the last section: 



1. OUTPUT.DAT 



computed cexact = 

r= .100000000D-01 

r= .200000000D-01 

r= .300000000D-01 

r= .400000000D-01 

r= .500000000D-01 

r= .600000000D-01 

r= .700000000D-01 

r= .800000000D-01 

r= .900000000D-01 

r= . 100000000D+00 

r= .110000000D+00 

r= . 120000000D+00 

r= . 130000000D+00 

r= . 140000000D+00 

r= . 150000000D+00 

r= . 160000000D+00 

r= . 170000000D+00 

r= . 180000000D+00 

r= . 190000000D+00 

r= .200000000D+00 

r= .210000000D+00 

r= .220000000D+00 

r= .230000000D+00 

r= . 240000000D+00 

r= .250000000D+00 

r= .260000000D+00 

r= .270000000D+00 

r= .280000000D+00 



.57142857D+00 
central charge= 
central charge= 
central charge= 
central charge= 
central charge= 
central charge= 
central charge= 
central charge= 
central charge= 
central charge= 
central charge= 
central charge= 
central charge= 
central charge= 
central charge= 
central charge= 
central charge= 
central charge= 
central charge= 
central charge= 
central charge= 
central charge= 
central charge= 
central charge= 
central charge= 
central charge= 
central charge= 
central charge= 



.571403656E+00 
.571329513E+00 
.571206952E+00 
.571036718E+00 
.570819520E+00 
.570556042E+00 
.570246948E+00 
.569892886E+00 
.569494492E+00 
.569052390E+00 
.568567192E+00 
.568039505E+00 
.567469925E+00 
.566859042E+00 
.566207440E+00 
.565515695E+00 
.564784379E+00 
.564014059E+00 
.563205295E+00 
.562358645E+00 
.561474660E+00 
.560553889E+00 
.559596875E+00 
.558604159E+00 
.557576277E+00 
.556513761E+00 
.555417143E+00 
.554286946E+00 



r= .290000000D+00 central charge= .553123695E+00 

r= .300000000D+00 central charge= .551927909E+00 
error in extrapolation .3365E-09 
estimated exponent .285714287D+01 
theoretical exponent .285714286D+01 
estimated dimension of the corresponding operator 
for a unitary theory: DELTA= .285714D+00 
for a non-unitary theory: DELTA= -.428571D+00 

fitted f_i 
f( 1)= .9643967331341316D-01 
f( 2)=-.1538311769518447D-02 
f( 3)= .6222166295705919D-04 
f( 4)=-.3197939259001748D-05 

chi-square value of the fitting= .5164E-29 

total cpu time (sees) .113E+06 



7 Conclusions 

We presented a multi-grid scheme for the resolution of the thermodynamic Bethe ansatz 
equations. The TBA is a means to describe the finite temperature effects of relativistic 
factorized scattering theories. Our program is specifically designed for theories having 
a scalar S-matrix. These theories exhibit a unique form, and the only input needed in 
order to carry out the TBA are the locations of the poles and zeros of the single S'-matrix 
elements. 

The program calculates the central charge and in the ultra-violet limit the dimension 
of the perturbing field and the coefficients of the perturbation expansion. These are 
the most crucial tests in verifying a conjectured S'-matrix. It should not be difficult 
for the user to add subroutines calculating other physical quantities, as for example for 
magnetic systems the moments of the total magnetization, or the convergence-region of 
the perturbation series in ([It]) |fj. [|. 



In order to get sensible results for the physical quantities one needs to resolve the 
integral equations with the highest possible accuracy. This unfortunately renders the 
calculation extremely time consuming. Therefore the use of an efficient Multi-Grid al- 
gorithm gives the possibility to reach high accuracy in the computation together with a 
sensible reduction of the CPU time, in confrontation with standard iterative techniques. 
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